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Abstract 

We present a numerical study of a two dimensional granular medium 
consisting of hard inelastic disks. The evolution of the medium through- 
out a cooling process is monitored. Two different types of instabilities 
(shearing and clustering instability) are found to develop in the sys- 
tem. The development of these instabilities is shown to be in qual- 
itative and quantitative agreement with the predictions of linearized 
hydrodynamic theory. 

1 Introduction 



Granular fluids are dense media composed of elementary elements of macro- 
scopic size, undergoing collisions in which their macroscopic energy is not 
conserved. In the last decade, the flow of this granular fluids has received a 
great attention from the physics community, both because these media offer 
a "simple" example of dissipative systems and because of their numerous 
industrial applications. Two factors make the behaviour of granular fluids 
very different from that of molecular fluids. Firstly, the macroscopic size of 
the particles implies that external fields (boundaries or gravity) have a much 
stronger effect on granular fluids. Secondly, the energy of the granular fluid 
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is not a conserved variable, since the heat dissipated in collisions can be 
considered as lost as far as the flow is concerned. These two effects are often 
difficult to disentangle in experiments, and also in the numerical simulations 
that aim at a realistic modelling of these experiments O. A different type 
of numerical simulations was initiated by several groups [0, |4], |5[ M. In 
this approach external fields are ignored, and only the dissipation is taken 
into account. This dissipation is moreover modelled in a very simplistic way, 
by describing the particles as monodisperse rigid hard spheres undergoing 
inelastic collisions. The dissipation is entirely specified by the restitution co- 
efficient r, where 1 — r is the fraction of the kinetic energy lost in a collision 
(in the center of mass frame) . 

The aim of these simulations is not to model actual experiments involving 
granular fluids. Instead, the models are used to assess the difference in be- 
haviour induced by the dissipation between a granular fluid and its " atomic" 
counterpart (r = 1, the hard sphere fluid.) In particular, the simulations 
can be used to investigate the validity of the hydrodynamic equations fre- 
quently used to describe granular flow in more realistic situations [0, g|. As 
for atomic fluids, they also provide a direct way of measuring the equation of 
state and transport coefficients that enter these equations. These transport 
coefficients can then be compared to those obtained using "granular kinetic 
theory" g. 

A particularly simple and instructive situation that can easily be studied 
in numerical simulations is the 2-dimensional "cooling problem" first stud- 
ied in H], H|. In their simulations, these authors start from an equilibrium 
configuration of an elastic (r = 1) hard disc fluid. The behaviour of this fluid 
after introducing a nonzero restitution coefficient is followed using the stan- 
dard molecular dynamics method for hard bodies [1C], with a collision rule 
that takes the dissipation into account. The kinetic temperature (average 
kinetic energy per particle) decreases due to the inelasticity of the collisions, 
so that the fluid cools down as time increases. It was shown in [||, ^] that, 
depending on the system size, on the restitution coefficient and density, this 
cooling follows different routes. In the simplest case, (small systems or small 
dissipation) the fluid remains homogeneous at all temperatures. In larger 
systems or for larger dissipations, either the velocity field or the density 
field in the fluid develop instabilities and become inhomogeneous. It was 
also found by the same authors that the occurence of such instabilities is in 
qualitative agreement with the predictions of a linear stability analysis of 
the hydrodynamic equations for granular fluids J?]]. Finally, it was discov- 
ered in H |5| that in some cases, the cooling ends at a finite time due to a 
singularity in the system dynamics, which was shown to correspond to an 
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infinite number of collisions within a finite time. This singularity, observed 
both in 1 and 2 dimensions, was described as an "inelastic collapse" of the 
system. 

In this paper, a detailed and quantitative analysis of the instability of 
homogeneous cooling of granular fluids is attempted. Our aim is to compare 
quantitatively the predictions of granular kinetic theory and granular hy- 
drodynamics to the results of molecular dynamics simulations of the cooling 
of an inelastic hard disk fluid. The paper is organized as follows. The main 
predictions of granular hydrodynamics concerning the cooling problem are 
briefly recalled. Computational details concerning the simulation are given 
in section 3. The different regimes occuring during the cooling are analyzed 
in section 4, and compared with the theoretical predictions. Our main focus 
will be on the growth rate of the density instability, that can be computed 
by monitoring the structure factor of the system as a function of time. Fi- 
nally, the problem of the "inelastic collapse" is adressed in section 5, where 
a possible method for avoiding this singularity in the system dynamics is 
proposed. 



2 Hydrodynamic analysis of the cooling problem 

The hydrodynamic equations that have been proposed to describe granular 
flow [|7j are based on mass and momentum conservation, and are very similar 
to the usual Navier-Stokes equations. The only modification is the appear- 
ance of a new term in the energy (or temperature) equation, accounting 
for the loss of energy in the collisions. These equations can be compactly 
written in the form 



g = (1, 
"15? = -VP ^ (2) 
= -V.Q-(r-(pB)- 7 rl (3) 

where D/Dt is the hydrodynamic derivative, D the symmetrized velocity 
gradient tensor, P the stress tensor, Q the heat flux and 7 represented the 
rate of energy lost due to inelastic collisions. For a hard disk fluid, the 
equation of state and the expression of the various transport coefficients can 
be obtained from Jenkins and Richman kinetic theory S. These expressions 
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are recalled in Appendix A. The energy sink term, 7T 3 / 2 , has also been 
written in the form appropriate for hard disks. 7 in that case is a function 
of the density and the restitution coefficient, which at least in the low density 
limit must be proportionnal to p and (1 — r). This can be understood from 
the following reasoning: the kinetic energy loss per particle per unit time is 
proportional to the collision frequency (i.e. to pT 1 / 2 ) and to the energy loss 
per collision (1 — r)T. 

A trivial solution of the cooling problem formulated above corresponds to 
an homogeneously cooling fluid, with a uniform density, a vanishing velocity 
field, and a uniform temperature with an algebraic time decay 

T(i)=T (l + £) 2 (4) 

. Here to = 2po/ '{iqTq 1 ^ 2 ) sets the time scale for temperature decay in the 
fluid. The linear stability of this homogeneous solution has been investigated 
in references pj. For completness, the main steps of this analysis will 
be repeated here. The linearized equations that describe the evolution of a 
sinusoidal perturbation 



5p = 5p k exp (ik • r) 
5v = 5v k exp (zk • r) 
ST = STk exp (zk • r) 

around the homogeneous solution are 
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As for usual fluids, the transverse part of the velocity field completely decou 
pies from the longitudinal part, and decays with time as (1 + t/to) 



fc 2 T 1/2 t /po 
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The longitudinal part of the velocity field, the temperature and the density 
are coupled, and give rise to three modes that have an algebraic time de- 
pendance 

5 Pk = 6p k [l + t/t ]t (12) 
5v k = 5^1 + t/tof- 1 (13) 

(14) 



5T k = 5T k [l + t/to]*- 



. The exponents for the three modes are the three roots of the deter- 
minant of the following set of equations 
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A typical plot of the wavevector dependance of these three roots, together 
with the growth rate of the velocity perturbations, is shown in figure 1. It 
must be emphasized that the stability of velocity disturbances is determined 
by the comparison between the growth exponent of the disturbance with 
the value —1 that characterizes the decay of the thermal velocity. Hence a 
growth exponent larger than —1 for the transverse or longitudinal velocity 
fields is indicative of an instability of the macroscopic velocity. If the growth 
exponent of the longitudinal velocity field is larger than —1, a corresponding 
instability in the density field will follow from equation ^. 

This analysis yields to the prediction of three different possible be- 
haviours of the system, depending on the value of the parameters and on 
the system size, that introduces a lower wavevector cutoff. If this lower cut- 
off corresponds to the line C of figure 1, the homogeneous solution will be 
linearly stable. This regime will be described as the homogeneous kinetic 
regime. If the lower cutoff moves to the abscissa indicated by line B in 
figure 1, the transverse velocity field will become unstable while the system 
remains homogeneous. In this "shearing" regime, first observed in reference 
H|, a shearing flow will develop in the system. Finally, for a lower cutoff 
corresponding to abscissa A, an instability of the longitudinal velocity field 
and the corresponding instability in the density field will take place together 
with the shearing instability. In this "clustering" regime, the growth of den- 
sity disturbances will yield to the formation of dense clusters of particles, as 
first observed in reference |^2| . 

All three situations have already been observed in numerical simulations 
of the cooling in two dimensional granular fluids. The aim of the next sec- 
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tions will be to attempt a quantitative analysis of the behaviour of a cooling 
granular fluid, and to compare the results to the predictions summarized 
above. 

3 Computational details 

The model simulated in this work is in all respects similar to that studied in 
p|. The system is made up of iV hard inelastic disks of diameter d, in a 
square cell of size L with periodic boundary conditions. The cell size L sets 
the lower cutoff in wavevector space, k m i n = 2ir/L. A standard cell-linked 
Molecular Dynamics algorithm for hard bodies |l(| is used. In a first step, 
the system is equilibrated with a coefficient r equal to unity. At time t = 0, 
inelasticity is switched on and cooling starts, with an initial temperature 
To. The restitution coefficient enters through a simple modification of the 
standard collision rule between hard disks, the velocities of the two disks 
after a collision being given by 

ui' = u 1 --(l + r)[n-(u 1 -u 2 )]n (15) 

u 2 / = u 2 + i(l + r)[n-(u 1 -u 2 )]n (16) 

where the primes denote the quantities after collision and n is a unit vector 
along the centers line from particle 1 towards particle 2. The natural units in 
this problem are the particle mass m and diameter, and the thermal energy 
at t = 0, i.e. To. The corresponding time unit is r = (m/T) l / 2 a. The 
state of the system is defined by three dimensionless numbers, which are 
the reduced size L/a (or equivalently the reduced cutoff k* min = k m i n a), the 
reduced density p* = a 2 N/L 2 , and the restitution coefficient r. 

The state of the fluid during the cooling was monitored by a systematic 
computation of coarse-grained (hydrodynamic) density and velocity fields. 
The coarse graining is obtaining here from a division of the system into 100 
square subcells. Besides, statistical quantities characterizing the state of 
the system have also been systematically computed. These quantities are 
the momenta of the velocity distribution of individual particles, the pair 
correlation function g(r) for interparticle distance, and the structure factor 

S(k) = ^p k p_ k (17) 

This structure factor can be computed for all wavevectors compatible with 
the periodic boundary condition, of the form (n x ,n y )k m i n . As the system 
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is not in a stationary state, these quantities are time dependant. A large 
enough system is thus necessary to obtain reasonable statistics without time 
averaging. The values of N investigated in this work vary from N = 1600 
to N = 10000. 

4 results 

4.1 Kinetic regime 

According to the analysis of section [l], the kinetic regime corresponding 
to a stable homogeneous cooling will be observed (at a given density and 
restitution coefficient) for small enough systems. Such a situation allows 
a clear testing of some of the hypothesis of the kinetic theory description 
of the granular fluid. In particular, the pair correlation function and ve- 
locity distribution can be compared to that of an elastic hard disk fluid 
throughout the cooling process. The temperature decay can be monitored 
and compared to the theoretical prediction (equation ||) , and the decay time 
to (or equivalently the coefficient *y(p)) compared to the prediction of kinetic 
theory. 

The pair correlation of an homogeneously cooling granular fluid after 
the temperature has dropped by a factor of 10 is shown in figure 2. This 
comparison shows that the local structure of the cooling granular medium 
(which determines its equation of state) remains essentially identical to that 
of an equilibrium fluid. The study of the velocity distribution function shows 
that this distribution remains maxwellian throughout the cooling. 

This similarity between the structure and velocity distribution of the 
granular fluid and the usual hard disk fluid suggests that the kinetic theory 
of Jenkins || is applicable. This expectation is borne out by the study 
of the time dependance of the fluid temperature. As shown in figure 3, 
the temperature decay is perfectly described by equation |I[ The density 
dependance of the decay time to is compared in figure 4 to the prediction 
of kinetic theory (see appendix B). The agreement is extremely good, and 
suggests that all the transport coefficients appearing in the hydrodynamic 
equations can be estimated using this kinetic theory. 

4.2 Shearing regime 

If the restitution coefficient r decreases or if the size of the system increases, 
the hydrodynamic theory predicts a regime in which transverse fluctuations 
of the velocity field are unstable. This regime is indeed observed in the sim- 
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ulations, as shown in figure 5. A shear flow that corresponds to the smallest 
wavevector compatible with the periodic boundary conditions develops in 
the system. In this regime, the total kinetic energy of the system (which 
in that case is not the temperature, since the system has developped an or- 
dered flow pattern) appreciably deviates from equation [|, as shown in figure 
6. 



4.3 Clustered regime 

For even larger systems or smaller restitution coefficients, the cooling gran- 
ular fluid becomes inhomogeneous, as shown in figure 7. this spontaneous 
formation of density inhomogeneities (or clusters) was first observed in the 
simulations of the cooling problem by Goldhirsch and Zanetti and Young 
and McNamara 0, Two different explanations have been put forward 
to explain this cluster formation. The first one, found in 0, is to consider 
this cluster formation as a secondary instability of the shearing regime, due 
to the developpment of temperature and pressure gradients in the shear- 
ing regime. The second possible explanation is that cluster formation is 
directly related to the linear instability of the density modes predicted by 
hydrodynamic theory. 

In order to characterize quantitatively this clustering regime, the struc- 
ture factor S(k, t) of the system has been computed as a function of time 
and wavevector. The corresponding data is shown in figure 8. The growth 
of the density inhomogeneities results in the appearance of a low wavevector 
peak in the structure factor, that rapidly increases with time. According to 
hydrodynamics, the time dependance of S(k, t) should be algebraic, i.e. 

S(k,t) = 5(k,0)(l + -J (18) 



so that the ratio 

ln(S(k,t)) -ln( S(k,0)) 



2£(k) (19) 



should be independent of time. This ratio is plotted in figure 9 as a function 
of wavevector for different times. 2£(fc) seems to be reasonably indepen- 
dent of time, and its low wavevector value appears to be consistent with 
the prediction of linearized hydrodynamics. Hence the density instability 
can be interpreted as resulting from a linear instability of the homogeneous 
solution of the hydrodynamic equations. Note that it was recently observed 
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by McNamara and Young that the "clustering" fluid eventually develops 
for long times into an ordered flow pattern of the "shearing" type. This is 
also consistent with hydrodynamics, since the growth rate of the transverse 
velocity modes is positive. The description of the formation of this shearing 
flow in an inhomogeneous system, however, is beyond the possibilities of 
linearized hydrodynamics. 

5 Inelastic collapse and how to avoid it 

The inelastic collapse singularity was first observed by || |||] in simulations 
of unidimensional inelastic system. This collapse can be described as the 
appearance of an infinite number of correlated collisions between a few par- 
ticles, taking place in a finite time. The same phenomenon was observed in 
two dimensions by ||. It was shown that in that case the correlated colli- 
sions take place between a small number of essentially aligned particles, so 
that the unidimensional situation is practically reproduced. 

In order to avoid this inelastic collapse, a slightly modified collision rule 
between the particles can be introduced. At each collision, the relative 
velocity of the two particles is first computed according to the usual rule 
(equations and |l^) , then rotated by a small (less than 5 degrees) random 
angle. This can be justified by invoking the unavoidable roughness of actual 
solid particles, conservation of angular momentum being (virtually) ensured 
by a transfer to the internal degrees of freedom of the particles. As to 
inelastic collapse, the aim of this modified collision rule is to hinder the 
formation of correlated particle lines that cause this singularity. Indeed, 
inelastic collapse was not observed in the simulations where this "random" 
collision rule was used, while under the same conditions a system following 
the "deterministic" collision rule always underwent inelastic collapse (figure 
10). Hence inelastic collapse appears to be a pathology related to the use of 
purely specular collision rule between particles, rather than a characteristic 
of inelastic fluids. 

6 Conclusion and perspectives 

The main objective of this work was to assess the validity of the hydrody- 
namic description of granular fluids originally proposed by Jjj], and of the 
kinetic theory calculation of the associated transport coefficients. The study 
of the particularly simple "cooling fluid" case and of the associated insta- 
bilities provides an ideal benchmark for this description. The comparison 
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between numerical simulations and theoretical predictions in this simple case 
shows that the theory is quantitatively accurate. A similar conclusion was 
also reached in a recent study by McNamara and Young fll3| , who showed 
that the transitions between the different cooling regimes were correctly 
predicted by the theory. 

The description of the inelastic collapse phenomenon observed by Mc- 
Namara and Young is obviously beyond the possibilities of kinetic theory or 
hydrodynamics. It was shown that this phenmenon can easily be avoided 
by introducing a small amount of randomness in the collisions between par- 
ticles, similar to what would be caused by the natural roughness of granular 
particles. 

Obviously, a correct description of granular fluid cannot be achieved 
without a knowledge of the boundary conditions that must be used for the 
hydrodynamic equations. These conditions, and in particular those that cor- 
respond to the very important case of vibrating solid walls, are not known. 
Their determination, through the quantitative comparison of numerical sim- 
ulation and theory, will be the subject of future work. 
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A Expressions for the transport coefficients and 
the equation of state 



The Navier-Stokes like equations describing a granular fluid are: 

pV- v 

V ■ T 

V • Q - tr (Plf) - 7T§ 

D/Dt is the hydrodynamic derivative, D the symmetrized velocity gra- 
dient tensor, P the stress tensor, Q the heat flux and 7 represents the rate 
of energy lost due to inelastic collisions. The definition for these quantities 
is: 



Dp 

Dt 
Dv 

7 Dt 
DT 

} Di 



D 9 , „. 

Dt = ^ + (V - V) 

Dij = + 

13 2 \dxj dxi) 

T = Ph T-2p(l5-^(V-v)T 

Q = - K VT 

The various transport coefficients and equation of state are 



Ph = p'{u) P T 

. 1 
p = p (ujTzpa 

k = U{v)T^pa 

1 = 1 (")- 
a 

where p', p! , k', 7' are functions only of the solid fraction v = pj ' p s . 



p'(v) = 



2v + s* 
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H\v) = 
K » = 

i{y) = 
where s* is defined as 

< \ {l -" r 
S * {U) = (1 - 7./16) 

B Enskog expansion 

In the case of a hard core fluid, a semi-empirical modification of Boltzmann 
equation introduced by Enskog, widens the range of applicabity of the ki- 
netic approach to higher densities [jlj] . The Enskog approximation accounts 
for the finite size of the disks in the collision term of the Boltzmann equation. 
When two particles collide, their centers are separated by the diameter of 
the disks a. The collision term of Boltzmann equation should thus be mul- 
tiplied by the probability of finding two particles separated by a which is 
proportional to the pair correlation function evaluated at a. This correction 
will have an influence on all the transport coefficients. The validity of this 
approach was checked for the cooling rate 7 in the kinetic regime. 

Figure 2 shows successive snapshots of the pair correlation function in the 
kinetic regime. This function is essentially the same as in a hard core fluid 
at thermodynamic equilibrium even though the temperature has dropped by 
a factor 10 between the first and the last snapshot. Hence g{a) is assumed 
to be given by the usual virial expression 

% = 1 + 2^) 

Introducing the equation of state of an 2d hard disks fluid |jTq| , the 
Enskog corrected cooling rate to becomes: 

1 

^0 ^0 Boltzmann 7 C 

I 1 o \ (l-^) 2 

\^={\-r)vT^) (1-7Z//16) 
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The values of the cooling rate found in the simulation are compared in 
figure 4 with this prediction. 
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figure 1: The decay rate of the velocity field disturbances computed by 
solving equation (4) for p = 0.2, r = 0.9 is shown in this figure as a function 
of wavevector. The dashed line indicates the decay rate of the transverse 
modes. The black and grey lines correspond respectively to the real and 
imaginary parts of the decay rate of the three longitudinal modes. 

figure 2: Pair correlation functions g(r) computed for the initial (in 
black shifted up by 0.5) and final (in grey) configurations in an homoge- 
neously cooling system (N = 1600, d = 0.8, r = 0.98). The temperature 
drops by a factor of ten without changes in the pair correlation function. 

figure 3: Evolution of the square root of the inverse temperature versus 
t in an homogeneously cooling system (r = 0.99, d = 0.1, N = 1600). The 
solid line corresponds to the hydrodynamic prediction in the kinetic regime. 

figure 4: Enskog corrected value of the decay time calculated as a 
function of density for r = 0.98, compared to the values obtained in the 
simulations (N = 1600). 

figure 5: Velocity field in the shearing regime after the granular medium 
has spontaneously developped a flow pattern corresponding to the lowest 
wave vector compatible with the boundary conditions. (N=1600, d=0.1, 
r=0.92) 

Ifigure 6: square root of the inverse of the temperature versus t in 
the shearing regime. The solid line extrapolates towards the first moments 
of the run. There is a substantial deviation from this kinetic regime fit. 
(N = 1600, d = 0.1, r = 0.92) 

figure 7 Final configuration (141 collisions per particle) of a simulation 
in the cluster regime. (N = 1600, d = 0.25, r = 0.6). 

figure 8: Evolution of the structure factor during the clusters formation 
(N = 1600, d = 0.5, r = 0.4). The curves from the bottom are separated by 
10 collisions per particle. Note the large increase of the structure factor in 
the long wavelength limit (k — > 0). 

figure 9: Growth exponent of the density field disturbance, obtained 
from the simulation using equation (19) (N=10000, d=0.5, r=0.9) The dif- 
ferent symbols correspond to different times. The prediction of the hydro- 
dynamic description of the instability is the solid line. 

figure 10 A system with TV" = 1600, d = 0.25, r = 0.25 obeying the 
specular collision rule collapses after 3.77 collisions per particle. The grey 
particles are those involved in the last two hundred collisions. The aligned 
particles represent more than 99 % of the grey particles. Under the same 
conditions, a system obeying the modified collision rule does not undergo 
collapse after 125 collisions per particle. 
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growth rate of the velocity disturbance 




This figure "fig5b.gif" is available in "gif" format from: 



http://arXiv.org/ps/cond-mat/9609184vl 
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